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Abstract. We review what has been learned recently using iV-body 
simulations about the evolution of globular clusters. While simulations 
of star clusters have become more realistic, and now include the evolu- 
tion of single and binary stars, the prospect of reaching large enough N 
is still a distant one. Nevertheless more restricted kinds of simulations 
CN \ have recently brought valuable progress for certain problems of current 

observational interest, including the origin and structure of tidal tails 
of globular clusters. In addition, such simulations have forced us to re- 
l^r) ' think some basic aspects of stellar dynamics, including, in particular, the 

■^j- \ process of escape. Finally we turn to faster, approximate methods for 

■ studying star cluster dynamics, where the role of iV-body simulations is 

one of calibration. 

7— I ! 

O 

O |- 1. Introduction 

2 ■ The title of this review unfortunately disguises the fact that it is about globular 

clusters, in the traditional sense. Actually, "a million" is somewhat on the large 
side in this context. Though the total number of stars in a globular cluster is 
not known, it may be estimated roughly by converting the luminosity (Harris 
1996) to the mass, using the method of Mandushev, Staneva, & Spasova (1991), 
and adopting a mean stellar mass of 0.34Mq , which corresponds approximately 
to a power-law mass function with slope —1 (Kroupa 2001) between turnoff and 
the hydrogen burning limit. The result is (Fig.l) that the median value is about 
3 x 10 5 , and roughly 10% of all clusters have more than one million stars. In this 
review we highlight some of the ways in which iV-body simulations are being 
used in studying the dynamical theory of globular star clusters. 

2. Comprehensive Calculations^ 

No simulation of a globular cluster will ever be complete, but rapid strides are 
being made to extend the traditional techniques of stellar dynamics so as to 
include stellar evolution. As a result there is a growing body of simulations 
with the following typical characteristics: (i) 1 star is represented by 1 particle 



At Tokyo several people expressed some distaste for the phrase "kitchen-sink calculations", 
which I used in my talk. 



1 



2 



Douglas C. Heggie 




100 1000 10000 100000 1e+06 1e+07 

N 



Figure 1. Cumulative distribution of approximate number of stars 
in the clusters of the galactic globular cluster system, derived as in the 
text. 

("star-by-star calculation"); (ii) forces include star-star gravity, a galactic tide 
(for a cluster on either a circular or elliptic galactic orbit), and disk shocking; 
(hi) treatment of stellar collisions; (iv) inclusion of stellar properties (luminos- 
ity, colour, metallicity) and stellar evolution; and (v) inclusion of (primordial) 
binaries and binary evolution. 

There is a strong synergy between the development of such codes and the 
development of special-purpose hardware (see the papers by Aarseth, Makino 
and the "starlab" group in this volume.) Table 1 lists some recent comprehensive 
simulations, where fun is the initial fraction of binaries. The scientific value of 
such work is well illustrated by Hurley et al (2001), who quantified the role of 
dynamical influences on the population of blue stragglers in an old open cluster. 
Such a result has been a long standing goal of this programme of research. 



Table 1. Some comprehensive calculations 



Author 


Object 


N(0) 


fbin 


Hurley et al (2001) 


M67 


15000 


50% 


Portegies Zwart et al (2001a) 


Pleiades, etc 


3000 


50% 


Portegies Zwart et al (2001b) 


Arches, etc 


12000 


0% 


Kroupa et al (2001)* 


"Pleiades" 


10000 


100% 



*No Roche Lobe mass transfer, etc 



Despite such successes, it is disappointing that not one of these is a sim- 
ulation of a globular cluster, in the traditional sense of the term. Even with 
GRAPE-4 it was possible to study a model with N = 32768 well past core col- 
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lapse (Makino 1996), which includes at least the smallest 10% of the existing 
galactic globular clusters. 

There are two reasons why simulation of any globular cluster is much harder 
than is suggested by the particle number alone. One is the requirement for a 
substantial population of primordial binaries. Though it has been argued that 
the computational load associated with binaries becomes negligible when N is 
sufficiently large (Makino & Hut 1990), this asymptotic limit is not yet in sight. 
The second reason is that the present particle number in a globular cluster may 
be a poor guide to what is required. Hurley et al (2001) estimated that even 
M67 had 7 times as many stars at birth as at present, and were therefore unable 
to simulate its early evolution. The present day mass function of the globular 
cluster NGC 6712 can be understood if it has lost about 99% of its original mass 
(Takahashi & Portegies Zwart 2000), and since its current mass is estimated to 
be of order 1O 5 M (Pryor & Meylan 1993), its likely original particle number 
puts it far beyond reach. 

One conclusion from this discussion is that, if a comprehensive simulation 
of a globular cluster is carried out within the next few years, it will tell us only 
about the evolution of the few clusters (if there are any) which are of low mass 
now and always have been. Study of the dynamical evolution of globular clusters 
cannot wait until such simulations can cope with all clusters. 

One of the aims of comprehensive simulations, as we have seen, is the study 
of the interaction between stellar and dynamical evolution, and the effects of 
collisions. While the dynamics of point masses can be computed with preci- 
sion, the effect of collisions is very uncertain. As Sills et al (2001) make clear, 
the effects of rotation on a merger remnant are dramatic. Since this is likely 
to be the largest source of uncertainty in studying the effects of collisions in 
globular clusters, it may be perfectly acceptable to adopt a model in which the 
point-mass dynamics is treated approximately, such as a Monte Carlo model 
(Sec. 5). These are much faster than direct iV-body methods, and can reach 
much larger "particle" numbers, ensuring that rare (but observable) species are 
not neglected. 



3. Restricted Simulations 

Considerable progress has been made in recent years in understanding important 
aspects of the dynamics of globular clusters without attempting a comprehensive 
simulation. The example we shall concentrate on here is the simulation of tidal 
tails. It is a theoretical problem that has been stimulated by observations (Grill- 
mair et al 1995, 1996; Kharchenko, Scholz, k, Lehmann 1997; Leon, Meylan, & 
Combes 2000; Testa et al 2000; Odenkirchen et al 2001). For this purpose one 
need include only (i) point mass dynamics, (ii) relaxation, (hi) the galactic tide, 
and (iv) the galactic orbit. Binaries and stellar evolution can be ignored. This 
is a pure stellar dynamics problem, and several sets of simulations have been 
carried out: Moore (1996); Combes, Leon & Meylan (1999), Murali & Dubinski 
(1999), Johnston, Sigurdsson & Hernquist (1999), Dehnen (pers. comm., and 
Fig.2), and others. The list would be much longer if one added comparable work 
on satellite galaxies. 
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Figure 2. Simulation of Pal 5 (Dehnen, pers. comm., reproduced by 
kind permission of the author). The orbit chosen is consistent with 
the present position, radial velocity and proper motion, and used the 
galactic potential of model 3 in Dehnen & Binney (1998). The initial 
condition is a King model with N = 10 4 . Upper panel: position on 
the sky (7,6); lower panel: radial velocity against longitude. Note the 
clumps in the tail in the upper panel; they resemble those in Pal 5 
itself (Odenkirchen et al 2001, where they are attributed to recent disk 
shocks). 
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What is not always clear from simulations is the mechanism by which stars 
escape from a cluster into the tidal tail. The most obvious possibilities are time- 
dependent tides (disk and bulge shocking) and relaxation. But such simulations 
usually cover only one or two Gyr of evolution, and it is possible that the choice 
of initial conditions plays a role, just as in the case of circular orbits (Fukushigc 
& Heggie 2000). The other hard thing to get right is relaxation, though the 
simulations of Combes et al (1999) show that this can be done. 

The role of these various mechanisms can be considerably clarified by a 
judicious choice of initial conditions and parameters. Baumgardt (this volume) 
begins at apogalacticon, with a cluster whose radius corresponds to the Roche 
radius at perigalacticon. The rationale for this is not the assumption that, over 
the course of many galactic orbits, the radius of a cluster is set by its Roche 
radius at perigalacticon; such an assumption would be in contradiction with 
observational evidence (Ninkovic 1985, Allen & Martos 1988, Meziane &: Colin 
1996; Brosche, Odenkirchen & Geffert 1999). Rather, if the initial radius of the 
cluster is either much greater or much smaller than the perigalactic tidal radius, 
then it is clear that initial conditions have much to do with the subsequent mass 
loss. 

Baumgardt's finding from his iV-body simulations is that, if one estimates 
the lifetime on the basis of mass lost in the first few orbits, as is often done, 
one finds that the lifetime increases with N roughly as the relaxation time, 
until sufficiently large N is reached, and then the lifetime is approximately 
independent of N. This is what would be expected if relaxation dominates 
for small iV and tidal heating dominates for large TV. If, however, one waits 
for as many galactic orbits as are required to lose a certain fraction (such as 
25%) of the mass, one finds that the mass loss is dominated by relaxation, for 
all N up to the largest values he was able to study. Though this N is still much 
smaller than the number of stars in the larger globular clusters, Baumgardt's 
iV-body models used particles of equal mass, and their relaxation time actually 
corresponds to the relaxation time of a much larger cluster with unequal masses. 

One draws two conclusions from this study. One is that the initial conditions 
matter; it is only after the lapse of a sufficient number of galactic orbits that the 
mass loss rate reaches its asymptotic size. Second, the escape rate is dominated 
by relaxation, and not by bulge shocking, even on a very eccentric galactic orbit. 
It remains to be seen, however, whether disk shocking is ever important. 

4. New Stellar Dynamics 

While it is clear that there are quantitative questions which can only be answered 
from iV-body simulations, has anything qualitatively new come out of them? 
Indeed, yes, but the list is not a long one. Most recently, we have learned that 
the lifetime of clusters does not scale with N in the expected manner, and the 
present section summarises this discovery. 

The aim of a collaborative experiment carried out in 1997 (Heggie et al 1998; 
Heggie, this vol.) was to compare results of various kinds of codes (and not only 
iV-body codes) in studying the dynamical evolution of a specific object. Since 
the number of particles was of order 250000, it was necessary to scale iV-body 
simulations by N. When this was done by the relaxation time, it was found 
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that the lifetime of the system depended on N, though its internal evolution 
(e.g. time to core collapse) was fairly insensitive to N. 

Fukushige suggested (pers. comm.) that the problem lay in the time taken 
for stars to escape, once they had enough energy to do so. If, however, it was 
assumed that the escape time scale is proportional to the crossing time, t cr , then 
it seemed hard to understand how the problem persisted for the large values of 
iV (up to 64K) for which results were available. The lifetime of a cluster should 
vary nearly with the relaxation time, except for N up to a few hundred, as in 
the escape theory of King (1959). Fukushige & Heggie (2000), however, found 
that the time taken to escape is a sensitive function of the energy, E, of a star: 
t esc oc (E — E esc )~ 2 , where E esc is the energy needed to escape (i.e. the potential 
at the Lagrange point, where the potential well of the cluster opens out into the 
rest of the galaxy). 

The next question was to find the appropriate energy to use in this formula. 
The answer is provided in Baumgardt (2001a). (His argument is more detailed 
than the following, but leads to the same scalings.) Unescaped stars with en- 
ergies just above E esc are relaxing during the time they take to escape, and so 
their excess energy is E — E esc oc (tesc/^rte) 1 , since relaxation is a diffusion 
process. (Here, t r i x is the relaxation time.) By solving these two relations we 
find the following conclusions: (i) t esc ~ (tutrix) 1 ^ 2 much longer than the cross- 
ing time, and quite different from previous expectations; this helps to explain 
why the effects of the escape time persist for much larger N than expected; (ii) 

E — E esc cx t~iJ^; this is also the fraction of stars inside the cluster with energy 
above E esc ; this fraction is still a few percent, even for systems of the size of a 
typical globular cluster; and (hi) from the previous two conclusions, the time to 

lose some fraction (say, half) of the mass of the cluster is tu oc t^. 

The last is the most far-reaching conclusion. Though the theory applies only 
for circular galactic orbits, Baumgardt has shown empirically (these proceedings) 
that the same conclusion appears to hold for clusters on elliptic galactic orbits. 
And yet Vesperini & Heggie (1997) studied the evolution of the mass function 
of clusters assuming that t\j oc t r i x , which is now known to be wrong. Indeed a 
great deal of recent modelling of globular cluster systems (e.g. Capuzzo-Dolcetta 
& Tesseri 1997, Gnedin & Ostriker 1997, Murali & Weinberg 1997, and other 
work of authors already quoted) may require revision. In effect, the new scaling 
of tu means that the lifetime of a globular cluster through dynamical processes 
is smaller (for large N) than one would have thought. Baumgardt (2001b) shows 
that one consequence of this is that the initial conditions of the cluster system 
are more completely erased. 

5. Calibration of Faster Methods 

5.1. Faster methods for star cluster dynamics 

Besides iV-body simulations, there are several other methods for studying the 
dynamical evolution of globular clusters (Fig.3). There are scaling methods, in 
which the cluster is characterised by nothing more than its mass and length 
scale (for example), Fokker-Planck methods, gas methods, and various hybrids. 
Even tree methods may function well for some important problems (McMillan & 
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Figure 3. Models of the dynamics of dense stellar systems. We 
do not include all hybrids, or methods used in other areas of stellar 
dynamics. 



Aarseth 1993, Arabadjis & Richstone 1998, Miocchi & Capuzzo-Dolcetta 2001). 
All these methods are faster and involve more approximations, and iV-body 
simulations are required in order to calibrate them. 

We have already seen (Sec. 2) that such methods may be quite appropriate 
for problems where the accuracy of the stellar dynamics is not the main uncer- 
tainty. Fast methods are also appropriate for the following kinds of problems: 

1. Evolution of cluster systems: For references see Sec. 4. For this purpose it 
must be possible to study the dynamical evolution of many globular cluster 
models with the full range of realistic N. At present this is possible only 
with fast approximate methods. 

2. Modelling individual objects: Fokker-Planck and (occasionally) fluid meth- 
ods have been used to construct evolutionary models for comparison with 
observations of specific clusters, including M3 (Angeletti, Dolcetta & Gian- 
none 1980), M71 (Drukier, Fahlman & Richer 1992), NGC 6624 (Grabhorn 
et al 1992), NGC 6397 (Drukier 1993, 1995), M15 (Grabhorn et al 1992, 
Phinney 1993, Dull et al 1997) and M30 (Howell, Guhathakurta & Tan 
2000). This requires extensive trial and error, with variations of the initial 
conditions of the model, and fast methods are essential. One motivation 
for such work has been the need to select suitable observational fields in 
clusters such that the local mass function requires least correction to the 
global mass function. 

3. Predicting initial conditions: This is the theorists' counterpart of the pre- 
vious point. If one wants to construct an iV-body model of a given object, 
considerable trial and error is needed in order to find appropriate initial 
conditions. Faster methods should be able to give a rapid but approximate 
answer, which can then be refined with a minimal number of iV-body sim- 
ulations. 
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5.2. Recent examples 

Scaling models As part of a study of the dynamical evolution of the mass func- 
tion of a star cluster, Vesperini & Heggie (1997) gave a simple formula for the 
evolution of the mass, M(t), i.e. M(t) = M(0) - AM s . e . - 0.828M(0)t/ F cw , 
where AM s . e . is the contribution from stellar evolution (which also depends 
on the mass function), and Few is the "family" parameter of Chernoff & 
Weinberg (1990). For a given mass function, it is proportional to the relax- 
ation time for average conditions inside the tidal radius, and is defined by 
Few = (M/M Q )(R g /kpc)(220kms~ l /v g )j 3 / 2 /InN, where M is the (original) 
mass of the cluster, moving on an orbit of radius R g at speed v g . 7 depends 
on the distribution of mass in the galaxy, is unity for a point-mass galaxy, and 
(3/2) - 1 / 3 for an isothermal potential. 

The above formula for M{t) has been used by Vesperini (1998, 2000, 2001) 
in studies of the evolution of cluster systems. Its relevance to the present review 
is that it was based on the simplest scaling considerations, but the numerical 
coefficient was obtained by examination of A-body simulations. In other words, 
A-body results were used to calibrate a simpler theoretical model. 

We now briefly consider the lifetime of a cluster. In the absence of stellar 
evolution, the above formula gives a cluster lifetime tyH = i 7 bw/0-828 in Myr. 
The coefficient depends weakly on the structure of the initial model. A similar 
formula for M{t) has been given by Portegies Zwart et al (2001b), which leads 
to a lifetime tpzMMH = 0.29t rX f. Here t rx t is another measure of the relaxation 
time within the tidal radius, defined by t rx t = 0.138M 1//2 r^ 2 /(mG 1 / 2 log 10 A) 
(cf. Spitzer 1987, eq.(2-63), except that Spitzer uses the natural logarithm in 
the Coulomb factor log A). In this formula, fh is the mean stellar mass. 

For the galactic mass law used by Portegies Zwart et al, 7 = (3/1. 8) 1 / 3 , and 
so it is easy to show, if we equate the two forms of the Coulomb logarithm, that 
the ratio of the lifetimes is tyH /t pzmmh — O.34M /m. The dependence on fh 
arises because Few /trxt also depends on fh. If it is assumed that the lifetime is 
proportional to t rxt (independent of the mass function, cf. de la Fuente Marcos 
1995), then it is necessary to modify tyH to tyu = {Few /0.828)(m\/_ff /fn), 
where fhyu is the value in the simulations by Vesperini & Heggie (1997). In 
fact for most of their work, which used a mass function f{m) oc m~ 2 - 5 for 
O.1M < m < 15M we have m V H — O.28M . The effect of this change 
(which is also necessary in the formula for M(t)) is to bring the two independent 
estimates of cluster lifetime into fair agreement. It must be borne in mind, 
however, that the assumption that the lifetime is proportional to the relaxation 
time is incorrect (Sec.4), and that fh is time-dependent. 

Despite the fact that lifetime (in the absence of stellar evolution) does not 
depend much on the initial structure (if the cluster initially fills its Roche lobe), it 
would be useful to have an equally simple way of estimating the time-dependence 
of the concentration parameter. 

Fokker-Planck models The comparison between A-body and Fokker-Planck 
models (computed with finite differences) has a considerable history, but recent 
work appears to have brought the two methods into reasonable agreement, with 
one caveat, discussed below. The central issue is the escape criterion, as shown 
by Takahashi & Portegies Zwart (1998, 2000). These papers show how a free 
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parameter in the Fokker-Planck specification can be determined by calibration 
using the results of iV-body models. Incidentally, the second of these is also 
an extensive survey of globular cluster models, which essentially supersedes the 
comparable extensive Fokker-Planck work of Chernoff & Weinberg (1990). 

As Takahashi & Portegies Zwart themselves point out, the caveat in this 
is that the iV-body models used in this survey adopt a tidal cutoff. Since it is 
known that the lifetime varies with N in a different way when a tidal field is used 
(Sec. 4), it remains to be seen how well the agreement between Fokker-Planck 
and iV-body models survives in the more realistic case. 

Monte Carlo models Though these can be regarded also as Fokker-Planck mod- 
els, they offer some advantages over finite-difference models, and have received 
increasing attention in recent years. While part of this effort involves calibra- 
tion with iV-body models, the real purpose of this little subsection is to say 
something about the present state of play. 

The largest models to date are those of Freitag &, Benz (2001), who present 
results of Monte Carlo simulations with up to 2 x 10 6 shells. At such levels 
it becomes possible to represent many globular clusters with one shell per star. 
These models, which are motivated by problems of galactic nuclei, do not include 
binaries, but have been extended to include stellar collisions (Fig. 4). For models 
including the effects of dynamically formed binaries the record is held, at time of 
writing, by Giersz (this volume). Reaching well past core collapse with N = 10 6 
shells, as he does, is a notable landmark. Primordial binaries are excluded 
from these models, but a 10% population of primordial binaries are present in 
computations with up to 3 x 10 5 shells by Rasio's group (Rasio, Fregeau & Joshi, 
2001; Fig. 5). Soon they will compute few-body interactions "on the fly", instead 
of using cross sections. All these developments give great promise for the Monte 
Carlo method in the coming years. 

6. Conclusions 

This rather selective review has concentrated on what iV-body models are telling 
us about the dynamics of globular star clusters. Comprehensive simulations, 
which include details of evolution of single and binary stars, are at an exciting 
stage, but are still limited to rich open clusters. It seems unlikely that com- 
prehensive simulations of globular clusters will be possible in the near future. 
Much progress can be made on some interesting kinds of dynamical problems 
about globular clusters (e.g. tidal tails) using more restricted kinds of simula- 
tion, which are perfectly feasible now. One of the dynamical lessons that has 
been learned in recent years is that it is vital to implement a correct treatment 
of tidal effects, even in the case of a steady tide (i.e. a circular galactic or- 
bit). Monte Carlo methods now have the dynamical aspects of globular clusters 
within reach. When a detailed treatment of stellar evolution has been incorpo- 
rated (at the level which has already been successful in the best iV-body codes), 
these codes may well prove optimal for the comprehensive study of globular star 
clusters, until the time when real iV-body models are fast enough. 

Acknowledgments. I thank the organisers of the meeting for their sup- 
port, and I am grateful for comments from H. Baumgardt, S.P.F. Portegies Zwart 
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Figure 4. Comparison of the collision rate between stars of mass M\ 
and M.2 in a static Plummer model (left) and in a Monte Carlo sim- 
ulation (right) with 512000 shells, but without stellar evolution (from 
Freitag 2000, with kind permission). The collision rate among more 
massive stars is greatly enhanced by mass segregation. 
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Figure 5. Simulation of an isolated system with N = 3 x 10 5 ini- 
tially, including 10% primordial binaries (F. Rasio, by kind permis- 
sion). Radii plotted against time are, from the bottom: core radius, 
half-mass radius of binaries, half-mass radius of all stars. Core col- 
lapse is avoided for about 70 initial half-mass relaxation times, when 
gravothermal oscillations begin. 
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and E. Vesperini on scientific matters. I also thank W. Dehnen, M. Giersz, M. 
Freitag and F. Rasio for the provision of unpublished results. The picture of M15 
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